************************************************************************
* 				MAIN REGRESSION TABS w/ FE and ME				       *
************************************************************************

use "data/output/iraqmerged.dta", clear

global treatsect i.treat##i.shia
global indvars c.age i.gender i.coledu
global disvars c.pctshia c.sqsect c.clanperthou
global regints i.regid i.regid#i.treat
global treatindcovints i.treat#c.age i.treat#i.gender i.treat#i.coledu

ebalance treat shia age gender coledu pctshia sqsect clanperthou if kurd==0 , targets(1) 


**********BALANCED FE*********

drop ident
clonevar ident=natident
logit ident $treatsect $indvars $disvars $treatindcovints $regints if kurd==0 [pw=_webal]
est sto m1nf
estadd local regdum "\checkmark"
estadd local regint "\checkmark"
estadd local treatintcovs "\checkmark"

drop ident
clonevar ident=arabident
logit ident $treatsect $indvars $disvars $treatindcovints $regints if kurd==0 [pw=_webal]
est sto m1af
estadd local regdum "\checkmark"
estadd local regint "\checkmark"
estadd local treatintcovs "\checkmark"

drop ident
clonevar ident=islident
logit ident $treatsect $indvars $disvars $treatindcovints $regints if kurd==0 [pw=_webal]
est sto m1if
estadd local regdum "\checkmark"
estadd local regint "\checkmark"
estadd local treatintcovs "\checkmark"

**********BALANCED ME*********

drop ident
clonevar ident=natident
melogit ident $treatsect $indvars $disvars $treatindcovints if kurd==0 [pw=_webal] || regid: ||  disid: treat shia
est sto m1nm
estadd local regslope "\checkmark"
estadd local treatintcovs "\checkmark"

drop ident
clonevar ident=arabident
melogit ident $treatsect $indvars $disvars $treatindcovints if kurd==0 [pw=_webal] || regid: ||  disid: treat shia
est sto m1am
estadd local regslope "\checkmark"
estadd local treatintcovs "\checkmark"

drop ident
clonevar ident=islident
melogit ident $treatsect $indvars $disvars $treatindcovints if kurd==0 [pw=_webal] || regid: ||  disid:
est sto m1im
//estadd local regslope "\checkmark"
estadd local treatintcovs "\checkmark"

***TABLE 1***

esttab m1nf m1af m1if m1nm m1am m1im using tables/iraqmlbFEME1.tex, se starlev(† .1 * .05 ** .01 *** .001) nonotes ///
 nobaselevels eqlabels("" "") varlabels(_cons "Constant" pctshia "\% Shi'i")  collabels("\textit{Coef.}" "\textit{SE}") ///
 mtitle("Nat. ident." "Arab. ident." "Islamic ident."  "Nat. ident." "Arab ident." "Islamic ident.") ///
 wide label replace cells("b(star fmt(3)) se(fmt(3))") scalars("regdum Province fixed effects" "regint Province*treatment" "treatintcovs Treatment*individual controls" "regslope Treatment/sect random slopes" ) ///
 drop(*.regid *.regid#*.treat *.treat#*.age *.treat#*.gender *.treat#*.coledu) noobs ///
addnotes("{\textit{† p$<$.1, * p$<$.05, ** p$<$.01, *** p$<$.001}}")


**********MARGINS PLOTS*********

***FIGURE 8 TOP PANEL***

local modlist = "m1nf m1af m1nm m1am"

foreach model of local modlist {

est resto `model'
margins, over(shia) at(treat=0) post noestimcheck
est sto bm
est resto `model'
margins, over(shia) at(treat=1) post noestimcheck
est sto am
	
coefplot (bm) (am), ///
	vertical recast(bar) barwidth(.3) level(68) ///
	ciopts(recast(rcap) color(gs10)) citop ///
	addplot(scatter @b @at, ms(i) mlabel(@b) mlabpos(2) mlabcolor(black)) format(%9.2f) ///
	xlab(1 "Sunni" 2 "Shi'i'") ///
	ytitle("") ylab(0(.2)1, format(%2.1f)) scale(.8) ///
	legend(order(1 "Before Mosul" 3 "After Mosul") size(small)) ///
	scheme(s1mono) aspectratio(1)  name(pr`model', replace)
gr_edit plotregion1.plot1.style.editstyle area(linestyle(color("%0"))) editcopy
gr_edit plotregion1.plot3.style.editstyle area(linestyle(color("%0"))) editcopy
gr_edit plotregion1.style.editstyle boxstyle(linestyle(color(none))) editcopy
gr_edit plotregion1.style.editstyle boxstyle(linestyle(width(none))) editcopy
}


grc1leg prm1nf prm1nm prm1af prm1am, ycommon rows(1) name(g1, replace)
gr_edit plotregion1.graph1.title.text.Arrpush National identity FE
gr_edit plotregion1.graph2.title.text.Arrpush National identity ME
gr_edit plotregion1.graph3.title.text.Arrpush Arab identity FE
gr_edit plotregion1.graph4.title.text.Arrpush Arab identity ME
gr_edit plotregion1.graph1.title.style.editstyle size(medium) editcopy
gr_edit plotregion1.graph2.title.style.editstyle size(medium) editcopy
gr_edit plotregion1.graph3.title.style.editstyle size(medium) editcopy
gr_edit plotregion1.graph4.title.style.editstyle size(medium) editcopy


**********BALANCED ME ADD INTERCEPT + TREAT INTS *********

global treatcovintsall i.treat#c.age i.treat#i.gender i.treat#i.coledu i.treat#c.pctshia i.treat#c.sqsect i.treat#c.clanperthou
 
drop ident
clonevar ident=natident
melogit ident $treatsect $indvars $disvars $treatcovintsall if kurd==0 [pw=_webal] || regid: ||  disid: treat shia
est sto m1nmri
estadd local regslope "\checkmark"
estadd local treatintcovsall "\checkmark"

drop ident
clonevar ident=arabident
melogit ident $treatsect $indvars $disvars $treatcovintsall if kurd==0 [pw=_webal] || regid: || disid: treat shia
est sto m1amri
estadd local regslope "\checkmark"
estadd local treatintcovsall "\checkmark"

drop ident
clonevar ident=islident
melogit ident $treatsect $indvars $disvars $treatcovintsall if kurd==0 [pw=_webal] || regid: || disid:
est sto m1imri
//estadd local regslope "\checkmark"
estadd local treatintcovsall "\checkmark"

**********BALANCED ME ADD COVS *********
 
global indvarsall c.age i.gender i.coledu i.wealthcatr i.voted i.hajjbin
global disvarsall c.pctshia c.sqsect c.clanperthou c.popdens i.urban c.sqmosdist
global rqualvars c.nmis c.intlsecs

lab var nmis "\#Missing responses"
lab var intlsecs "Completion time"

drop ident
clonevar ident=natident
melogit ident $treatsect $indvarsall $disvarsall $treatcovintsall $rqualvars if kurd==0 [pw=_webal] || regid: || disid:
est sto m1nmriac
estadd local treatintcovsall "\checkmark"

drop ident
clonevar ident=arabident
melogit ident $treatsect $indvarsall $disvarsall $treatcovintsall $rqualvars  if kurd==0 [pw=_webal] || regid: || disid:
est sto m1amriac
estadd local treatintcovsall "\checkmark"

drop ident
clonevar ident=islident
melogit ident $treatsect $indvarsall $disvarsall $treatcovintsall $rqualvars  if kurd==0 [pw=_webal] || regid: || disid:
est sto m1imriac
estadd local treatintcovsall "\checkmark"

***TABLE B.1***

esttab m1nmri m1amri m1imri m1nmriac m1amriac m1imriac using tables/iraqmlbMEac.tex, se starlev(† .1 * .05 ** .01 *** .001) nonotes ///
 nobaselevels eqlabels("" "") varlabels(_cons "Constant" pctshia "\% Shi'i")  collabels("\textit{Coef.}" "\textit{SE}") ///
 mtitle("Nat. ident." "Arab ident." "Islamic ident." "Nat. ident." "Arab ident." "Islamic ident.") ///
 wide label replace cells("b(star fmt(3)) se(fmt(3))") scalars("treatintcovsall Treatment*individual/district controls" "regslope Treatment/sect random slopes") ///
 drop(*.treat#*.age *.treat#*.gender *.treat#*.coledu *.treat#*.pctshia *.treat#*.sqsect *.treat#*.clanperthou) noobs ///
addnotes("{\textit{† p$<$.1, * p$<$.05, ** p$<$.01, *** p$<$.001}}")

**********MARGINS PLOTS*********

***FIGURE B.4***

local modlist = "m1nmri m1amri m1nmriac m1amriac"

foreach model of local modlist {

est resto `model'
margins, over(shia) at(treat=0) post noestimcheck
est sto bm
est resto `model'
margins, over(shia) at(treat=1) post noestimcheck
est sto am
	
coefplot (bm) (am), ///
	vertical recast(bar) barwidth(.3) level(68) ///
	ciopts(recast(rcap) color(gs10)) citop ///
	addplot(scatter @b @at, ms(i) mlabel(@b) mlabpos(2) mlabcolor(black)) format(%9.2f) ///
	xlab(1 "Sunni" 2 "Shi'i'") ///
	ytitle("") ylab(0(.2)1, format(%2.1f)) scale(.8) ///
	legend(order(1 "Before Mosul" 3 "After Mosul") size(small)) ///
	scheme(s1mono) aspectratio(1)  name(pr`model', replace)
gr_edit plotregion1.plot1.style.editstyle area(linestyle(color("%0"))) editcopy
gr_edit plotregion1.plot3.style.editstyle area(linestyle(color("%0"))) editcopy
gr_edit plotregion1.style.editstyle boxstyle(linestyle(color(none))) editcopy
gr_edit plotregion1.style.editstyle boxstyle(linestyle(width(none))) editcopy
}


grc1leg prm1nmri prm1nmriac prm1amri prm1amriac, ycommon rows(1) name(g3, replace)
gr_edit legend.style.editstyle boxstyle(linestyle(width(vthin))) editcopy
gr_edit legend.plotregion1.DragBy 26.5 1.4 editcopy
gr_edit plotregion1.graph1.title.text.Arrpush National identity ME 2
gr_edit plotregion1.graph2.title.text.Arrpush National identity ME 3
gr_edit plotregion1.graph3.title.text.Arrpush Arab identity ME 2
gr_edit plotregion1.graph4.title.text.Arrpush Arab identity ME 3
gr_edit plotregion1.graph1.title.style.editstyle size(medium) editcopy
gr_edit plotregion1.graph2.title.style.editstyle size(medium) editcopy
gr_edit plotregion1.graph3.title.style.editstyle size(medium) editcopy
gr_edit plotregion1.graph4.title.style.editstyle size(medium) editcopy

graph export "plots/prmarginsacs.pdf", replace

**********MARGINS PLOTS ALL*********

local modlist = "m1nf m1af m1nm m1am m1nmri m1amri m1nmriac m1amriac"

foreach model of local modlist {

//make separate predictions for shia and sunni at treat=0/1
est resto `model'
margins r.treat, over(r.shia) noestimcheck post
est sto mc`model'
}

coefplot (mcm1nf,label(FE)) (mcm1nm,label(ME)) (mcm1nmri,label(ME 2)) (mcm1nmriac,label(ME 3)), bylabel({bf:Nat. ident.}) || ///
mcm1af mcm1am mcm1amri mcm1amriac, bylabel({bf:Arab. ident.}) ||, ///
vertical yline(0, lcolor(black) lwidth(vvthin)) title(,size(tiny)) ytitle("Contrasts, Shi'i : Sunni", size(vsmall)) title("") ///
		levels(95) ylab(,nogrid labsize(vsmall)) xlab("",nogrid labsize(vsmall)) ciopts(lwidth(thin)) ///
		scheme(s1mono) ///
		subtitle(, size(small)) ///
		mlabposition(2) mlabsize(vsmall) ///
		mlabel format(%12.2f) mlabposition(3.5) mlabgap(*2) mlabsize(vsmall) ///
		legend(rows(1) size(small)) byopts(row(1)) title("Model specification", size(3.5)) aspectratio(.3) name(g2, replace)
gr_edit plotregion1.title[1].text = {}
gr_edit plotregion1.title[2].text = {}


gr combine g1 g2, rows(2)
gr_edit plotregion1.graph1.legend.style.editstyle boxstyle(linestyle(color(none))) editcopy
gr_edit plotregion1.graph2.legend.style.editstyle boxstyle(linestyle(color(none))) editcopy
graph export "plots/prmarginmcs.pdf", replace

**********ADDITIONAL COV P-STAT LOOPS*********
**********BALANCED NAT*********

tuples i.wealthcatr i.voted i.hajjbin c.popdens c.urban c.sqmosdist
di `ntuples'*1
// prepare an empty matrix in which to store the
// p-values
matrix pval = J(`ntuples', 1, .)
global base c.age i.gender i.coledu c.pctshia c.sqsect c.clanperthou
 
forval i = 1/`ntuples' {
    cap melogit natident $treatsect $indvars $disvars $treatindcovints `tuple`i'' if kurd==0 [pw=_webal]  || regid: || disid:
		if _rc==0 {
				margins treat#shia, contrast post
				scalar p = el(r(table),4,1)
				scalar vars = e(est_cmdline)
				matrix pval[`i++',1] = p
				di vars
				}
				else if _rc!=0 continue 
}

matlist pval
putexcel set "data/output/pvalsnb.xlsx", sheet("pval") replace
putexcel A1=matrix(pval)

preserve
import excel "data/output/pvalsnb.xlsx", sheet("pval") clear
rename A pvalue
histogram pvalue, freq bin(20) ///
color(gs6) lcolor(gs6) ///
xtitle("{it:p}-value", size(small )) ///
ytitle("Frequency", size(small)) ///
xlab(, labsize(small) format(%9.4f)) ///
ylab(, labsize(small)) ///
scheme(s1mono) ///
ylab(, nogrid) ///
title("Nat. ident.", size(small)) name(natb, replace) aspectratio(.5)
graph export "plots/tryoutmcsnbnew.png", replace
restore


**********BALANCED ARAB*********

tuples i.wealthcatr i.voted i.hajjbin c.popdens c.urban c.sqmosdist
di `ntuples'*1
// prepare an empty matrix in which to store the
// p-values
matrix pval = J(`ntuples', 1, .)
global base c.age i.gender i.coledu c.pctshia c.sqsect c.clanperthou
 
forval i = 1/`ntuples' {
    cap melogit arabident $treatsect $indvars $disvars $treatindcovints `tuple`i'' if kurd==0 [pw=_webal]  || regid: || disid:
		if _rc==0 {
				margins treat#shia, contrast post
				scalar p = el(r(table),4,1)
				scalar vars = e(est_cmdline)
				matrix pval[`i++',1] = p
				di vars
				}
				else if _rc!=0 continue 
}

matlist pval
putexcel set "data/output/pvalsab.xlsx", sheet("pval") replace
putexcel A1=matrix(pval)

preserve
import excel "data/output/pvalsab.xlsx", sheet("pval") clear
rename A pvalue
histogram pvalue, freq bin(20) ///
color(gs6) lcolor(gs6) ///
xtitle("{it:p}-value", size(small )) ///
ytitle("Frequency", size(small)) ///
xlab(, labsize(small) format(%9.4f)) ///
ylab(, labsize(small)) ///
scheme(s1mono) ///
ylab(, nogrid) ///
title("Arab ident.", size(small)) name(arb, replace) aspectratio(.5)
graph export "plots/tryoutmcsabnew.png", replace
restore


gr combine natb arb, rows(1)

graph export "plots/tryoutmcsall.pdf", replace
